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ABSTRACT 

We present the results of a numerical study examining the effect of H2-based star formation (SF) 
model on the rest-frame UV luminosity function and star formation rate function (SERF) of z > 
6 galaxies, and the implications for reionization. Using cosmological hydrodynamical simulations 
outfitted with an H2-SF model, we find good agreement with our previous results (non-H2 SF model) 
and observations at M^v < —18. However at M^v > —18, we find that the LF deviates from both our 
previous work and current observational extrapolations, producing significantly fewer low luminosity 
galaxies and exhibiting additional turnover at the faint end. We constrain the redshift evolution of 
this turnover point using a modified Schechter function that includes additional terms to quantify the 
turnover magnitude (M*^) and subsequent slope (/3). We find that M*^ evolves from M*^ = —17.33 
(at 2: = 8) to —15.38 (z = 6), while /? becomes shallower by A/3 — 0.22 during the same epoch. This 
occurs in an Afuv range which will be observable by James Webb Space Telescope. By integrating the 
SERF, we determine that even though H2-SF model significantly reduces the number density of low 
luminosity galaxies at Afuv > —18, it does not suppress the total SER density enough to affect the 
capability of SF to maintain reionization. 

Subject headings: cosmology: theory — galaxies: evolution — galaxies: formation — methods: numer- 
ical 



1. INTRODUCTION 

Characterizing the properties of galaxy population at 
high redshift is one of the most important goals of mod- 
ern cosmological study, as star-forming galaxies are the 
primary candidate for the sources of ionizing photons 
that reionized the universe. In recent observational (e.g., 
Bouwens et al. 2012a; Finkelstein et al. 2012; McLure 
et al. 2012; Schenker et al. 2012; Oesch et al. 2012) and 
theoretical (e.g., Trenti et al. 2010; Salvaterra et al. 2011; 
Jaacks et al. 2012a; Kuhlen & Faucher-Giguere 2012) 
works, characteristics of the z > 6 UV LF have been mea- 
sured and predicted with increasing accuracy. They both 
suggest that the faint-end slope of UV luminosity func- 
tion (UV LF) is quite steep, and that the ionizing pho- 
ton budget is dominated by the low-mass, star-forming 
galaxies below the current detection limit. 

In particular, some cosmological hydrodynamic simu- 
lations have found extremely steep faint-end slopes of 
a < —2.0, which could be the result of too efficient star 
formation (SF) in their subgrid SF model at high-z. Pre- 
vious SF models (e.g., Springel & Hernquist 2003; Schaye 
& Dalla Vecchia 2008; Choi & Nagamine 2010) relied on 
the statistical representation of locally observed Kenni- 
cutt (1998) relationship, which is dependent on the to- 
tal gas column density. Recent observations (Kennicutt 
et al. 2007; Leroy et al. 2008; Bigiel et al. 2008) sug- 
gest however that SF correlates strongly with molecu- 
lar hydrogen (H2). As the H2 gas is less abundant at 
higher redshifts due to lower metallicities, the previous 
SF models that rely only on total gas density may be 
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overproducing stars (Krumholz et al. 2009; Kuhlen et al. 
2012), potentially affecting the number density of low- 
mass galaxies at high-z. 

In this work, we quantify the effect of H2-SF model 
on our simulated UV LF, star formation rate function 
(SERF), and total SER density (SERD). This paper is 
organized as follows: Section 2 briefly discusses the meth- 
ods and parameters of our simulations; Section 3 presents 
our results, and we summarize in Section 4. 

2. SIMULATIONS 

We use a modified version of the smoothed particle 
hydrodynamics (SPH) code GADGET-3 (originally de- 
scribed in Springel 2005). Our code includes radia- 
tive cooling by H, He, and metals (Choi & Nagamine 
2009), heating by a uniform UV background (UVB; 
Faucher-Giguere et al. 2009), the UVB self-shielding 
effect (Nagamine et al. 2010), the initial power spec- 
trum of Eisenstein & Hu (1999), supernova feedback, 
sub-resolution multiphase interstellar medium model 
(Springel & Hernquist 2003), and the Multicomponent 
Variable Velocity (MVV) wind model (Choi & Nagamine 
2011). Our current simulations do not include active 
galactic nucleus feedback. 

Simulations are setup with either 2 x 400^^, 2 x 500^^ or 
2 X 600'^ particles for both gas and dark matter. Multiple 
runs were made with comoving box sizes of 10/i~"'^Mpc, 
34/i^^Mpc and 100/i^^Mpc to cover a wide range of halo 
and galaxy masses. Throughout this work they will be re- 
ferred to as N400L10, N500L34 and N600L100 runs. See 
Table 1 for simulation details. The adopted cosmolog- 
ical parameters are consistent with the latest WMAP7 
data: fl^ = 0.26, Qa = 0.74, f^b = 0.044, h = 0.72, 
as = 0.80, Us = 0.96 (Komatsu et al. 2011). We note 
that the Fiducial runs used the Salpeter (1955) stellar 
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Table 1 

Simulation Parameters Used in this Paper. The parameter Np is the 
number of gas and dark matter particles; ttidm and the particle 

masses of dark matter and gas; e is the comoving gravitational softening 

length. 



Run 


Box Size 


Np 


muM 




e 




(/i-lMpc) 


(DM,Gas) 




(/i-iM©) 


(/i-lkpc) 


N400L10 


10.0 


400^ 


9.37x10^ 


1.91x10^ 


1.0 


N500L34 


33.75 


500^ 


1.84x10'' 


3.76x10^ 


2.70 


N600L100 


100.0 


600^ 


2.78x10* 


5.65x10^ 


4.30 



initial mass function (IMF), while the H2 runs used the 
Chabrier (2003) IMF for historical reasons of our work. 
Galaxies are identified and grouped based on the bary- 
onic density field (see Nagamine et al. 2004, for more 
details) . 

Since the estimation of H2 mass fraction is dependent 
upon metallicity, the details regarding our feedback and 
enrichment models are relevant. When SF takes place, 
metals arc also prodiiccd with an instantaneous yield of 
0.02, and thereafter tracked by the code based on a closed 
box model for each gas particle (i.e., no diffusion). Our 
MVV wind model is designed to account for both energy- 
driven and momentum-driven winds (Choi & Nagamine 
2011). Wind velocity is determined by 7Jwi„d = C^'csc. 

where v^sc = 130(SFR)i/3 (i±2)^/^ kms"!. We adopt 
the standard values of ( = 1.5 for high-density regions 

(momentum driven) and C — 1 for low-density regions 
(energy driven), chosen by Choi & Nagamine (2011). 
The mass loading factor is = (cro/a-gai)^ for the energy- 
driven case, and rj = CTo/ugai for the momentum-driven 
case, where ctq = SOOkms"-^ and cTgai — Wosc/2 is the 
velocity dispersion of a galaxy. For full detail and phys- 
ical justifications for this model, see Choi & Nagamine 
(2011). 

Our "Fiducial" runs use the "Pressure-SF model" 
(Schaye & Dalla Vecchia 2008; Choi & Nagamine 
2010). while the present work uses the H2-SF model of 
Krumholz et al. (2009), implemented by Thompson ct al. 
(2013, hereafter H2 runs). This equilibrium analytic 
model calculates the SFR based on the II2 mass density 
rather than the total cold gas density, and Krumholz & 
Gnedin (2011) have shown that it is in good agreement 
with more computationally expensive, non-equilibrium 
calculations by Gnedin et al. (2009). The details of the 
implementation and the basic results of this model have 
been presented by Thompson et al. (2013). 

In principle our implementation of the H2-SF model of 
Krumholz et al. (2009) must be similar to the previous 
work by Kuhlen et al. (2012) on the most basic level. 
The primary difference between the two work is in the 
class of code in which it was implemented, Enzo (AMR) 
versus GADGET (SPH). We will further discuss the ba- 
sic differences and potential effects in Sections 3.1 and 
4. 

3. RESULTS 

3.1. Modified Schechter Luminosity Function 

We combine the results of our three runs to create a 
composite LF, which covers a much wider dynamic range 
than is possible with a single cosmological run. In Fig- 
ure 1, we present our composite LF for z = 6,7,8 (red 



triangles, blue circles, green squares) for the H2 run, in 

comparison to the Schechter (1976) fits for our Fiducial 
runs (dashed red, blue, green lines; Jaacks et al. 2012a) 
with the Pressure-SF model. We also show the observed 
LF fit range (gray shade; Bouwens et al. 2011). A small, 
constant extinction Eb-v is required to fall within the 
observational range for both runs (Jaacks et al. 2012a), 
although the II2 runs at 2; = 7 & 8 require less extinction 
by AEb-v = 0.025 than the Fiducial runs, suggesting a 
trend of decreasing Eb~v with increasing redshift. 

The value of Eb-v is chosen to be consistent with the 
value used to match the observed rest-frame UV LF in 
our previous work (Jaacks et al. 2012a), and it is centered 
between the following two recent observations: Bouwens 
et al. (2012b) argued for little to no extinction at the 
faint end of the LF at z = 6, whereas Willott et al. (2012) 
found a best-fit value of Ay = 0.75, which corresponds 
to Eb-v ^ 0.19 assuming Ry = 4.05 (Calzetti et al. 
2000) at the bright end of the UV LF at 2; = 6. This 
moderate amount of extinction is also consistent with the 
estimates by Schaerer & de Barros (2010) and de Barros 
et al. (2012) who included nebular emission lines in their 
spectral energy distribution fits. Therefore the values 
of Eb-v chosen for this work are reasonably consistent 
with current observations. 

At Muv ^ —18, both Fiducial and H2 runs show excel- 
lent agreement with each other and observations. How- 
ever at Muv > — 18, the II2 run start to show a turn-over 
of the LF, which is not present in the Fiducial LFs. This 
flattening significantly reduces the number density of low 
luminosity objects in the H2 run, and it occurs beyond 
current observational limit of the Hubble Space Telescope 
(HST). As this population of low-luminosity galaxies is 
thought to be the critical contributor to the total ioniz- 
ing flux at these redshifts (Trenti et al. 2010; Salvaterra 
et al. 2011; Bouwens et al. 2012a; Finkelstein et al. 2012; 
Jaacks et al. 2012a), it is important to quantify this re- 
duction and its implications. 

To quantify the turnover point and flattening, we adopt 
a modified Schechter function (hereafter Schechter-|-): 



$(L) = 



exp - 



(1) 



where </>*, L* and a are the normalization, characteristic 
luminosity, and faint-end slope of the standard Schechter 
function. The additional parameter indicates the 
point at which the LF undergoes its second turn, and 
/3 is related to the power-law slope at the lowest lumi- 
nosities. Note that $(i) cx L"^'^ when L -C i*, and 
that both a and /3 take negative values. A similar func- 
tional form to Equation (1) was used by Loveday (1997) 
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Figure 1. Composite rest-frame UV LF at z = 6,7,8 for the H2 run, shown as red triangles, blue solid circles, and green triangles, 
respectively. The solid lines represent the Schechter+ fit (Equation(l)) at each redshift to these simulation points. We also show the fits 
to the Fiducial run results (dashed lines) obtained by (Jaacks et al. 2012a). The observed LF fit range (gray shade; Bouwens et al. 2011; 
Oesch et al. 2012) are also shown, which are extrapolated to Muv > —18 based on their published fits. Observed points from McLure 
et al. (2012) and Schenker et al. (2012) are shown by the open black circles and red crosses at ^ = 7 & 8, respectively. The vertical dashed 
line indicates the approximate HST observational limit {Muv ~ —18) at these redshifts and the vertical dash-dotted line represents the 
approximate JWST limit [Muv ~ —16). The z = 8 panel also includes 2; = 6 & 7 line (red & blue) to demonstrate evolution. We show the 
results of Enzo AMR simulation from Kuhlen et al. (2012, black open squares) in the 2 = 6 panel for direct comparison. 



to characterize the upturn in dwarf galaxy population at 
low redshift. 

Using the magnitude version of Equation (1), we per- 
form a least-square fit to the H2 runs by fixing the three 
standard Schechter parameters {(j>* , M*^,a) to those ob- 
tained for the Fiducial runs by Jaacks et al. (2012a), and 
vary the remaining two parameters The best- 

fit results can be seen as solid red, blue and green lines 
in Figure 1 and are quantified in Table 2. 

In the rightmost z = 8 panel, we have overplotted the 
fits to z = 6 & 7 results to demonstrate the evolution. 
One can clearly see that the number of galaxies increases 
from z = 8 to 2 = 6 as the dark matter structure develops 
and more SF takes place. 

One might note the somewhat large uncertainty in the 
estimation of (3. This is due to the fact that we are fix- 
ing the standard Schechter parameters {(j)* , M*^,a) to 
our Fiducial results and only fitting the new Schechter -|- 
parameters (Af*.^,, /3). This was done in order to make a 
valid comparison between the current results and our pre- 
vious results (Jaacks et al. 2012a). If we make a full five 
parameter fit, we noticed that the value of M*^ becomes 
unphysically brighter and a becomes quite steeper, due 
to the chosen Schechter-f functional form. As a result, 
the size of the error bar of /? becomes smaller, however 
there would be no meaningful physical comparison to 
be made between the three parameter Schechter values 
which was obtained in our previous work. Therefore we 
chose to fix the standard Schechter parameters as in our 
previous work, despite of large error bar on /3. 

In the leftmost z = 6 panel, we have over plotted the 
Enzo AMR simulation from Kuhlen et al. (2012, black 
open squares, corrected for dust extinction by the same 
amount as our simulation), which also used the same 
Krumholz et al. (2009) II2-SF model. At the brighter side 
of Muv < —18, we see good agreement between the two 
simulations. However, on the fainter side of A/uv > —18, 



Table 2 

Best-fit Schechter -|- LF Parameters at 2 = 6 — 8 (Top to Bottom) 
for the H2 Runs. The errors arc Icr. 



log(<^-) A/*, M^, 13 

(Fixed) (Fixed) (Fixed) 



2.151 


.24 
.15 


-3.46±;f, 


-21.15t 


53 
.53 


-15.38t;^60 


-i.sot 


47 
.47 


2. sol 


.28 
.18 


-3.74t;4i 


-20.82j; 


.61 
.56 


-16.25t}:33 


-1.291: 


.70 
.70 


2.5lt 


.27 
.17 


-4.30+;42 


-2i.ooj; 


.59 
.50 




-1.281 


.49 
.49 



the Enzo result drops sharply relative to our SPH re- 
sults. While the origin of this difference is not entirely 
clear, one possible explanation is that the AMR simu- 
lation is missing the early structure formation of low- 
mass dark matter halos relative to our SPH simulations 
(O'Shea et al. 2005). The simulation of Kuhlen et al. 
(2012) used a comoving box size of 12.5 Mpc, which is 
not very different from our N400L10 run. However, their 
root grid is only 256'^, which effectively corresponds to 
128^ particles of SPH/7V-body simulation (O'Shea et al. 
2005). Although their highest spatial resolution after re- 
cursive refinement might be higher than our simulation, 
their AMR simulation would miss a large number of low- 
mass halos at early times due to lack of force resolution 
at early times, hence significantly lower number of dwarf 
galaxies compared to our SPH simulations which used 
400'^ dark matter particles. This is a well-known issue of 
AMR cosmological simulation, and has been extensively 
discussed by O'Shea et al. (2005). If Kuhlen et al. (2012) 
wanted to resolve the low-mass halos down to the same 
halo mass as our N400L10 run, they had to start their 
simulation with ^ 800"^ root grid and correspondingly 
low dark matter particle mass. In their Figure 1, Kuhlen 
et al. (2012) show that their halo mass function starts to 
deviate from Sheth & Tormen (1999) at A/haio < lO^M©, 
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Figure 2. Redshift evolution of Scliechter+ parameters for botli 
our simulated LF (red squares) and simulated SFRF (blue dia- 
monds) at z = 6, 7, & 8, shown with la error bars. The SFRF /3 
value is offset by Az = 0.1 for readability. The horizontal dash- 
dotted line in the top panel represents the approximate JWST 
observable limit (Muv ~ —16). 



whereas our halo mass function in N400L10 run agrees 
very weU with Sheth & Tormen (1999) down to IO^Mq, 
as we showed in Figure 17 of Jaacks et al. (2012a). This 
corroborates our suspicion that the main difference in the 
two results is coming from the difference in the number 
of low-mass halos in the two simulations. We also dis- 
cuss the other possibility of difference in metal diffusion 
in Section 4. 

In Figure 2, we summarize the redshift evolution of 
best-fit Schechter-I- parameters. We see that M*^ be- 
comes dimmer and the power-law slope a — /3 of $(i) 
becomes shallower with decreasing redshift. These re- 
sults support the current paradigm of hierarchical struc- 
ture formation: as the smaller objects form first and later 
merge to form larger systems, the number of dim object 
decreases, and the LF becomes flatter at the faintest end. 
Since the predicted Af*^ occurs at Muv < — 16 at z = 7 
& 8, this feature should be observable by future missions 
such as James Webb Space Telescope {JWST). 

3.2. Star Formation Rate Function (SFRF) 

A recent work by Smit et al. (2012) presented the ob- 
served SFRF from dust-corrected UV LF by converting 
Ltjy into SFR using the Kennicutt (1998) conversion. 
The SFRD is then easily obtained by integrating the 
Schechter-like fit to the SFRF down to a certain lower 
limit, culminating in the commonly used Lilly-Madau 
diagram (Lilly et al. 1996; Madau et al. 1996). Compar- 
ing the SFRF of our simulated galaxies to the observed 
results is an excellent test of our SF model, since SFR is 
one of the most basic outputs of our simulations. It is 
also a more intrinsic comparison for us as it does not re- 
quire any assumptions on our part regarding the amount 
of extinction in our simulated galaxies as this correction 
is included in the observational estimates. 

Figure 3 shows the SFRF of simulated galaxies at 



z = 6, 7, 8 (red triangles, blue circles and green squares, 
respectively) . We find good agreement with the observed 
results (Smit et al. 2012, cyan diamonds), especially at 
z = 7. Deviation at low-SFR end of z = 6 is expected, 
given a similar deviation in the LF at the low-end shown 
in Jaacks et al. (2012a, b). This deviation could either 
indicate that our simulation is still overproducing stars 
at z = 6, or it could also be due to uncertainties in the 
faint-end observations (i.e., missed faint galaxies) and/or 
assumptions regarding extinction. 

The SFRF has a similar functional form as the LF, 
therefore we utilize the same Schechter-|- function to fit 
it: </>(SFR) = dn/dlog(SFR) = ln(10)SFR $(SFR), and 



</.(SFR) = ln(lO)0* ( ) exp 



VSFR* 



SFR 



SFR* 

/ SFR \ 
lySFRt J 



(2) 



where (f>* , SFR* and a are the usual Schechter parame- 
ters as in Equation (1); SFR* is the location of the second 
turnover, and /? is related to the low-SFR end power-law 
slope which is proportional to (l-t-a— /3) at SFR ^ SFR*. 
The best-fit results to all five parameters are shown in 
Figure 3 with solid lines. Similarly to the rest-frame UV 
LF, the right-most panel shows the increase in normal- 
ization and brightening of SFR* clearly from z = 8 to 
z = 6. 

3.3. SFRD and Reionization 

Having a continuous representation of the SFRF allows 
for easy integration and acquisition of the SFRD at each 
redshift. In Figure 4 we show the results of this integra- 
tion (red solid and green lines) plotted against observa- 
tional results by Bouwens et al. (2011, black squares) and 
Oesch et al. (2012, blue diamond). The red solid line rep- 
resents an integration down to log(SFRiim[M0yr~^]) = 
0.0, which is approximately consistent with the lowest 
SFR value of the observational data. The green solid 
crossed line is the SFRD value when integrated down to 
log(SFRiiin) ~ —1.0. This value is chosen to be consis- 
tent with a galaxy stellar mass Ms = 10'' via the Mg- 
SFR relationship found in our simulations (SFRoc ■°). 
The mass of Ms — IO^Mq also represents the minimum 
galaxy resolution of the N400L10 run (^100 star parti- 
cles). The solid observational data points were obtained 
assuming the Salpeter IMF, therefore we also show the 
same data points corrected for the Chabrier IMF with 
open symbols by applying a simple factor (see Section 
3.2.1 of Thompson et al. 2013). Our conclusions are not 
affected by the change in the IMF. 

Figure 4 demonstrates that, when integrated down to 
a minimum SFR (log(S'i^i?iim) = -0.10), the H2 run 
agrees very well with the Fiducial run and observational 
estimates. This is expected, since the H2-SF model did 
not affect the number density of objects within the ob- 
served range of Afuv < —18 (Figure 1). We note that the 
observed points were calculated by integrating the UV 
LF to a hmiting Muv ~ -17.74 Bouwens et al. (2012b), 
which corresponds to log(SFR) « —0.10 when converted 
via Kennicutt (1998) relation. 

The green crossed line, representing the integration 
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Figure 3. SFRF of simulated galaxies at 2 = 6, 7, & 8, shown as red triangles, blue circles and green squares, respectively. The observational 
estimates (Smit et al. 2012) are shown by filled cyan diamonds. Solid red, blue and green lines represent the best-fit Schechter-|- functions 
(Equation 2) to the simulation data. 
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Figure 4. SFRD obtained from the integration of SFRF for 
both H2 and Fiducial runs with SFR limits of log(SFR) = -0.10 
(solid/dashed red line) and log(SFR) = —1.0 (solid/dashed green 
crossed line). Observations from Bouwens et al. (2012b) are rep- 
resented by the black squares and recent z = 8 observations from 
Oesch et al. (2012) by a blue diamond. The open black squares 
and open blue diamond show the same observed points adjusted for 
a Chabrier IMF. The SFRD required to maintain IGM ionization 
(Equation 3) is shown by the shaded gray area for 7 < C/ fcsc < 50. 
The blue dot-dashed line represents an updated estimation of the 
critical SFRD by ShuU et al. (2012). 



limit of log(SFRiini) = -1.0, shows that the total SFRD 
in our simulations is significantly higher than what is 
currently observed. This is consistent with our previous 
findings (Choi & Nagamine 2012; Jaacks et al. 2012a). 
The reduction of SFRD at z > 6 is also consistent with 
findings by Krumholz et al. (2009); Gnedin & Kravtsov 
(2010), and Kuhlen et al. (2012), who show that H2- 
SF model reduces high-z SFRD due to metallicity effect 
(Thompson et al. 2013). However the degree of reduc- 
tion may still be different among different simulations 
and models. 

To determine whether or not our simulated galaxy pop- 



ulation is sufficient to maintain reionization, we utilize 
the theoretical prescription presented in Madau et al. 
(1999), which quantifies the minimum SFRD required 
to keep the intergalactic medium (IGM) ionized (shaded 
gray contour): 



« 2 X 10" 



/esc 



1 



10 



[M0yr-iMpc-3]. (3) 



This depends on redshift and the ratio of IGM clumping 
factor (C) and escape fraction (/esc) of ionizing photons 
from galaxies. Given that the exact values of both C and 
/esc are still debated and uncertain, we show a wide range 
of 7 < C//csc < 50 (grey shade in Figure 4) which are 
consistent with works by Iliev et al. (2006); Pawlik et al. 
(2009): Finlator et al. (2012); Kuhlen & Faucher-Giguere 
(2012). 

We also include an updated estimation of the critical 
SFRD by ShuU et al. (2012, blue dot-dashed line) which 
includes considerations for the Ly-continuum production 
rates and the temperature scaling of the recombination 
rate coefficient. Their calculation uses a fiducial value of 

C / fcsc — 15. 

Figure 4 reinforces our previous arguments (Jaacks 
et al. 2012a) that the abundant, low-luminosity galax- 
ies, which are currently beyond the detection threshold of 
HST, dominate the total SFRD at z > 6. When this pop- 
ulation is considered, there is sufficient amount of ioniz- 
ing photons available to maintain ionization by z = 6 — 7 
for a reasonable value of C//esc- 

4. CONCLUSIONS AND DISCUSSIONS 

Using GADGET-3 cosmological SPH simulations 
equipped with a II2-SF model, we examined UV LF, 
SFRF and the contribution of low-luminosity galaxies to 
the total SFRD at z > 6 . Our major conclusions are as 
follows. 

• We find that, at Muv < -18, the Il2-based SF 
model does not change the faint-end slope from 
our Fiducial runs with a < — 2.00atz > 6 
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(Jaacks ct al. 2012a). However at the fainter end of 
Muv > —18, we find a significant reduction in num- 
ber density of faint galaxies due to the metallicity 
effect of H2-SF model (Thompson et al. 2013). 

• To characterize the new LF's turnover magnitude 
(M*y) and subsequent slope (/3), we utihze the 
Schechter+ function (Equation 1), and find that 
Ml^ evolves with redshift from M*^ = -15.33 
{z = 6) to M*^ = -17.69 {z = 8), while the second 
power-law slope {1 + a — P) becomes steeper over 
the same redshift range with A/3 = 0.22 (Figure 2). 

• We present the SFRF of simulated galaxies (Fig- 
ure 3), and find good agreement with observations. 
We also find similar evolution in the Schechtcr+ pa- 
rameters of SFRF: as log(SFR*) becomes slightly 
lower with increasing redshift, and the second 
power- law {1 + a — P) becomes steeper over the 
same redshift with A/3 = 0.52 (Figure 2). 

• Although the H2-SF model reduces the number 

density of low-luminosity galaxies, there is still suf- 
ficient SFRD (i.e., ionizing photons) to maintain 
rcionization at 6 < z < 8 (Figure 4) for a wide pa- 
rameter range of clumping factor (C) and escape 
fraction (/esc)- 

When examining the results of numerical simulations, 
the resolution effect is always a concern. The N400L10 
run used in the present work is reaching the densities 
similar to that of giant molecular clouds (GMCs), but it 
is not resolving the formation and internal structure of 
GMCs adequately. The model of Krumholz et al. (2009) 
is meant to capture the general behavior of GMCs on av- 
erage, and Thompson et al. (2013) have shown that the 
H2-SF model has a number of advantages over the previ- 
ous SF models. Nevertheless, as discussed by Thompson 
ct al. (2013), the H2-SF model has some inherent resolu- 
tion dependency. Although we believe that the observed 
turnover in the UV LF and SFRF are physical in nature, 
the exact values of M*^ and log(SFR*) may vary with 
resolution. 

Another concern is the impact of UVB on the faint- 
end slope of LF. Earlier work using semi-analytic models 
(e.g.. Babul & Rces 1992; Bullock et al. 2000; Benson 
et al. 2003) and simulations (e.g., Thoul & Weinberg 
1996; Gnedin 2000; Hoeft et al. 2006; Okamoto et al. 
2008; Finlator et al. 2011) have shown that the UVB 
suppresses the formation of low-mass galaxies via pho- 
toionization of gas in shallow potential wells. In par- 
ticular, Okamoto et al. (2008) have shown that the gas 
fraction drops rapidly from cosmic fraction (r^b/f^m) at 
Mtot^lO^MQ to almost zero at A^tot^lO^Af©. It is there- 
fore plausible that the turnover seen in our LF and SFRF 
could be the result of photoevaporation of gas from low- 
mass halos and not a feature of the implemented H2-SF 
model. In order to verify that the turnover at M*^ and 
SFR* in the II2 runs arc not affected by the UVB, we 
ran two test simulations with and without UVB. As a 
result, we found that the turn-over still exists even in 
the no-UVB run, and only the normalization of UV LF 
decreased slightly due to photoionization. The overall 



shape of the LF and SFRF remained unchanged, i.e., 
the values of M*^ and SFR* did not change significantly. 
Therefore, while the UVB reduces the number of low- 
mass galaxies in our simulations, we confirmed that the 
turn-over at M*^ and SFR* is due to the II2-SF model. 

In Section 3.1 wc discussed the differences between the 
resuhs of Enzo AMR simulation and our GADGET SPH 
simulations. Another fundamental difference between 
the two codes is in metal mixing, which could have an ef- 
fect on these results. In our SPH code metals are locked 
to the individual SPH particles, thus limiting mixing and 
diffusion of metals. In contrast grid-based codes follow 
the fluid element in cells, which may allow for more nat- 
ural mixing and diffusion within and between the cells. 
Conversely in low resolution regions, this could lead to 
over-smoothing of metals. At the moment, it is unclear 
how much impact the metal mixing and diffusion would 
have on the results presented in this paper, unless we 
perform more rigorous comparison tests between the two 
codes. We intend to investigate the issue of metal mix- 
ing and diffusion in the future using both SPH and AMR 
codes. 
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